clear all
version 15
cd "_____"
set more off

local idval=67 // idval=7 for Beijing , idval=67 for Shanghai
forv yearval=2001/2010{
use "./data/pm10data.dta", clear
qui keep if id==`idval' // city by city
qui keep if year==`yearval' // year by year
scalar cff1=0.15
sum pmconc if pmconc<=float(cff1)
gen Fxnp=r(N)/_N  
gen cens=0

*****************
* CENSORED MLE *
*****************
* manipulation window W1 
scalar lftcff1=0.135   
scalar rtcff1=0.22
replace cens=1 if pmconc>float(lftcff1) & pmconc <=float(cff1)  
gb2lfit_mod pmconc, cdf(pmparacdf1)  censvar(cens)
local ahat=e(ba)
local bhat=e(bb)
local phat=e(bp)
local qhat=e(bq)
local atx=cff1
gen Fxp1_cff=ibeta(`phat',`qhat', (`atx'/`bhat')^`ahat'/(1+(`atx'/`bhat')^`ahat'))
local atx=rtcff1
gen Fxp1_rtcff=ibeta(`phat',`qhat', (`atx'/`bhat')^`ahat'/(1+(`atx'/`bhat')^`ahat'))
gen prop1=Fxnp-Fxp1_cff
gen inrange1=Fxp1_rtcff-Fxp1_cff
gen prop1_inrange=prop1/inrange1

* manipulation window W2 
scalar lftcff1=0.135   
scalar rtcff1=0.18
replace cens=1 if pmconc>float(lftcff1) & pmconc <=float(cff1)  
gb2lfit_mod pmconc, cdf(pmparacdf2)  censvar(cens)
local ahat=e(ba)
local bhat=e(bb)
local phat=e(bp)
local qhat=e(bq)
local atx=cff1
gen Fxp2_cff=ibeta(`phat',`qhat', (`atx'/`bhat')^`ahat'/(1+(`atx'/`bhat')^`ahat'))
local atx=rtcff1
gen Fxp2_rtcff=ibeta(`phat',`qhat', (`atx'/`bhat')^`ahat'/(1+(`atx'/`bhat')^`ahat'))
gen prop2=Fxnp-Fxp2_cff
gen inrange2=Fxp2_rtcff-Fxp2_cff
gen prop2_inrange=prop2/inrange2
*****************
* POLYNOMIAL *
*****************
qui gen xrnd=round(pmconc,0.005) // round x to 0.005. Chetty's method requires discrete variables
gen xdisc=0.005*_n-0.005 
replace xdisc=. if xdisc>1
gen xfreq=.
forvalues i=1/201 {
sum xrnd if xrnd==float((`i'-1)*0.005)
replace xfreq=r(N)/_N in `i'
}
gen xdisc2=xdisc^2
gen xdisc3=xdisc^3
gen xdisc4=xdisc^4
gen xdisc5=xdisc^5
gen xdisc6=xdisc^6
gen xdisc7=xdisc^7

* order 6, window 1
forvalues i=28/45{
local val=(`i'-1)*5
gen xdum`val'=(xrnd==float((`i'-1)*0.005))
}
reg xfreq xdisc xdisc2 xdisc3 xdisc4 xdisc5 xdisc6 xdum*
forvalues i=28/45{
local val=(`i'-1)*5
replace xdum`val'=0
}
predict xfreqhat
sum xfreqhat in 1/31
gen Fxppoly61=r(mean)*r(N)
gen pmpolycdf61=sum(xfreqhat)
drop xdum* xfreqhat

* order 5, window 1
forvalues i=28/45{
local val=(`i'-1)*5
gen xdum`val'=(xrnd==float((`i'-1)*0.005))
}
reg xfreq xdisc xdisc2 xdisc3 xdisc4 xdisc5 xdum*
forvalues i=28/45{
local val=(`i'-1)*5
replace xdum`val'=0
}
predict xfreqhat
sum xfreqhat in 1/31
gen Fxppoly51=r(mean)*r(N)
gen pmpolycdf51=sum(xfreqhat)
drop xdum* xfreqhat

* order 7, window 1
forvalues i=28/45{
local val=(`i'-1)*5
gen xdum`val'=(xrnd==float((`i'-1)*0.005))
}
reg xfreq xdisc xdisc2 xdisc3 xdisc4 xdisc5 xdisc6 xdisc7 xdum*
forvalues i=28/45{
local val=(`i'-1)*5
replace xdum`val'=0
}
predict xfreqhat
sum xfreqhat in 1/31
gen Fxppoly71=r(mean)*r(N)
gen pmpolycdf71=sum(xfreqhat)
drop xdum* xfreqhat

* order 6, window 2
forvalues i=28/37{
local val=(`i'-1)*5
gen xdum`val'=(xrnd==float((`i'-1)*0.005))
}
reg xfreq xdisc xdisc2 xdisc3 xdisc4 xdisc5 xdisc6 xdum*
forvalues i=28/37{
local val=(`i'-1)*5
replace xdum`val'=0
}
predict xfreqhat
sum xfreqhat in 1/31
gen Fxppoly62=r(mean)*r(N)
gen pmpolycdf62=sum(xfreqhat)
drop xdum* xfreqhat

* order 5, window 2
forvalues i=28/37{
local val=(`i'-1)*5
gen xdum`val'=(xrnd==float((`i'-1)*0.005))
}
reg xfreq xdisc xdisc2 xdisc3 xdisc4 xdisc5 xdum135-xdum180
forvalues i=28/37{
local val=(`i'-1)*5
replace xdum`val'=0
}
predict xfreqhat
sum xfreqhat in 1/31
gen Fxppoly52=r(mean)*r(N)
gen pmpolycdf52=sum(xfreqhat)
drop xdum* xfreqhat

* order 7, window 2
forvalues i=28/37{
local val=(`i'-1)*5
gen xdum`val'=(xrnd==float((`i'-1)*0.005))
}
reg xfreq xdisc xdisc2 xdisc3 xdisc4 xdisc5 xdisc6 xdisc7 xdum*
forvalues i=28/37{
local val=(`i'-1)*5
replace xdum`val'=0
}
predict xfreqhat
sum xfreqhat in 1/31
gen Fxppoly72=r(mean)*r(N)
gen pmpolycdf72=sum(xfreqhat)
drop xdum* xfreqhat
local cityname=city[1] 

cumul pmconc, gen(pmecdf) 
gen cutoff=pmconc*0+cff1

local cityname=city[1] 
line pmecdf pmconc if pmconc<=0.4, sort xlabel(0(0.1)0.4) ylabel(0(0.2)1) yscale(range(0 1.1)) lcolor(black) graphregion(color(white)) ///
|| line pmparacdf1 pmconc if pmconc<=0.4, lcolor(navy) sort ///
|| line pmparacdf2 pmconc if pmconc<=0.4, lcolor(navy) lpattern(dash) sort ///
|| line pmecdf cutoff, lcolor(black) lpattern(dot) title("`cityname', `yearval': Censored MLE, W1/W2")  ///
name(`cityname'`yearval'_MLE, replace) aspect(0.8) ///
legend(order(1 "ECDF" 2 "MLE-W1" 3 "MLE-W2") ///
rows(2) size(4) region(c(none)) bm(zero)) ytitle("") 

line pmecdf pmconc if pmconc<=0.4, sort xlabel(0(0.1)0.4) ylabel(0(0.2)1) yscale(range(0 1.1)) lcolor(black) graphregion(color(white)) ///
|| line pmpolycdf51 xdisc if xdisc<=0.4, lcolor(red) ///
|| line pmpolycdf61 xdisc if xdisc<=0.4, lcolor(orange) ///
|| line pmpolycdf71 xdisc if xdisc<=0.4, lcolor(purple) ///
|| line pmecdf cutoff, lcolor(black) lpattern(dot) title("`cityname', `yearval': Polynomial Fitting, W1")  ///
name(`cityname'`yearval'_POLYw1, replace) aspect(0.8) ///
legend(order(1 "ECDF" 2 "Poly5-W1" 3 "Poly6-W1" 4 "Poly7-W1") ///
rows(2) size(4) region(c(none)) bm(zero)) ytitle("") 

line pmecdf pmconc if pmconc<=0.4, sort xlabel(0(0.1)0.4) ylabel(0(0.2)1) yscale(range(0 1.1)) lcolor(black) graphregion(color(white)) ///
|| line pmpolycdf52 xdisc if xdisc<=0.4, lcolor(red) lpattern(dash) ///
|| line pmpolycdf62 xdisc if xdisc<=0.4, lcolor(orange) lpattern(dash) ///
|| line pmpolycdf72 xdisc if xdisc<=0.4, lcolor(purple) lpattern(dash) ///
|| line pmecdf cutoff, lcolor(black) lpattern(dot) title("`cityname', `yearval': Polynomial Fitting, W2")  ///
name(`cityname'`yearval'_POLYw2, replace) aspect(0.8) ///
legend(order(1 "ECDF" 2 "Poly5-W2" 3 "Poly6-W2" 4 "Poly7-W2" ) ///
rows(2) size(4) region(c(none)) bm(zero)) ytitle("") 
}

local cityname=city[1] 
graph combine  `cityname'2007_MLE  `cityname'2007_POLYw1 `cityname'2007_POLYw2, ///
rows(1) graphregion(color(white)) ysize(7) xsize(20)  
graph export "./outputgraphs/main/Figure3_`cityname'.pdf", replace


